home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Precision Software Appli…tions Silver Collection 4
/
Precision Software Applications Silver Collection Volume 4 (1993).iso
/
stats
/
chadyn.exe
/
ADDMAP.DOC
next >
Wrap
Text File
|
1988-10-15
|
9KB
|
280 lines
AN EXAMPLE OF HOW TO ADD A NEW MAP
The following two functions bob() and init_bob() were added
to YMAPS.C to add the map called by the name "rr".
bob() /* this routine defines the map */
{
int rand();/* a Microsoft random integer generator */
double temporary, c, s, angle;
if ( rand() < C2*32768.)/* rand() returns a random integer
between 0 and 32768, so C2 should be a number
between 0 and 1, so that the first process will be
carried out C2 of the time */
{
y[0] = 2. -y[0];
y[1] = -y[1];
}
else
{
angle = C1*twopi/360;
c = cos(angle);
s = sin(angle);
temporary = rho*(c*y[0] + s*y[1]);
y[1] = rho*(-s*y[0] + c*y[1]);
y[0] = temporary;
}
}
init_bob()/* this process initializes the constants and with the first
line tells the program the name of the routine that defines the
map; since C1, C2, and rho are used in bob(), htey must be
initiallized in this routine; if they are not, the program will
not let be accessed interactively and they will be given default
values, namely -9999. */
{
map = bob;/* this tells the program the name of the routine that
is evaluated; map is a pointer to routines and bob()
must be defined in the file prior to this line; */
rho = .65;
C1 = 135.;
C2 = .5;
/* the following 4 determine the x and y scale */
X_lower = -3.;
X_upper = 5.; /* x scale */
Y_lower = -3.;
Y_upper = 3.; /* y scale */
}
The following expressions were added to YMAPMENU.C into the
routine mapmenu(). The name "rr" is name of the routine that you will
use to access thses routines. Use lower case letters.
TESTMAP2 is used for expressions that will appear in two ways:
the menu and the header. The header appears in the Main Menu and the
Parameter Menu and the printout. The text (in quotes) includes the map
name so the user knows what to call the process. This first instruction
results in the following line:
RR at Random either Rotate by C1 degrees and mult by rho OR flip about 1,0
TESTMAP (no "2") is for material that appears only in the header.
Note that TESTMAP is used twice since two lines of header material
is to be used. THe result is the following header:
RR at Random either Rotate by C1 degrees and mult by rho OR flip about 1,0
The probability of a flip (x -> 2-x, y-> -y) is C2
Otherwise rotate C1 degrees about 0,0 and multiply x and y by rho (<1)
The source code for the program is:
TESTMAP2("rr") fputs(
"RR at Random either Rotate by C1 degrees and mult by rho OR flip about 1,0\n"
,output);
TESTMAP("rr")
fputs(
" The probability of a flip (x -> 2-x, y-> -y) is C2\n"
,output);
TESTMAP("rr")
fputs(
" Otherwise rotate C1 degrees about 0,0 and multiply x and y by rho (<1)\n"
,output);
The following was added to YMAPMENU.C into the routine init_map() so that
the map could be initiallized (using the initialization routine which I
called init_bob()); initialization is important in part because it tells
the program what the name of the function is (bob in this case);
TESTMAP("rr")
{
init_bob();
continue; /* causes exit from switch */
}
THE VECTOR y[] COMPUTING LYAPUNOV EXPONENTS
The vector y[] can have dimension up to EQUATIONS which is now 21.
The vector y includes the phase space variables and time for differential
equations and all the coordinates of the tangent vectors used in computing
Lyapunov exponents.
1: the dimension of the phase space (not including time); this is also
the number of coordinates of the Lyapunov vectors;
2: the number of Lyapunov vectors to be computed initially; this can be
changed while running;
3: lyapzero is the number of coordinates exclusive of the lyapunov vector
coordinates; the first coordinate of the first Lyapunov
vector is y[lyapzero];
The following are used in the routine initHenon in the file YHENON.C
to specify these numbers.
zeroth = 0; /* the first coordinate that is a phase space
coordinate; it is assumed that the "vec_dim" phase space
coordinates are together so that coordinate
zeroth + vec_dim
is the first coordinate following the phase space coordinates;
the non phase space coordinates may either correspond to
the time variable or the time modulo some number such as
twopi (which is defined) or may be some other variable of
interset such as the total sum of the values of some
coordinate;
vec_dim = 2; /* the dimension of the
Lyapunov vectors = phase space dim; this number is used in a
number of places; it is used in computing lyapunov exponents;
it tells what the dimension of the Lyapunov vector is;
it is similarly used in Newton method calculations; it is also
used when the distance between two vectors is needed; the
distance usually useds only the phase space coordinates */
zeroth = 0; /* the first coordinate that is a phase space
coordinate; it is assumed that the "vec_dim" phase space
coordinates are together so that coordinate
zeroth + vec_dim
is the first coordinate following the phase space coordinates;
the non phase space coordinates may either correspond to
the time variable or the time modulo some number such as
twopi (which is defined) or may be some other variable of
interset such as the total sum of the values of some
coordinate;
num_lyap = 0; /* This is set only if you give the program the
equations for computing Lyapunov exponents;
the number of Lyapunov numbers
to be computed <=vec_dim; if this is not set; you will not be
able to change its value later from within the program; if
you do not tell the program what the partial derivatives are
then you do not set num_lyap; the fact that it is not set to
any value means it takes the default value of -9999 */
lyapzero = 2; /* y[lyapzero] is the zeroth coord of the
zeroth lyapunov vector */
Computation of Lyapunov exponents is done in such a way that very
little is wasted when Lyapunov exponets are not wanted. The following
is a variant of the actual routine Henon() in file YHENON.C.
Henon()
{
int i;
int base;
double y_temp[6];
y_temp[0] = rho -y[0]*y[0] + C1*y[1];
y_temp[1] = y[0];
/* the following "for(...){...}" is for computing tangent vectors; */
for (i = 0; i < num_lyap; i++)/* if the number of Lyapunov
exponents to be computed num_lyap is currently 0,
then this part of the routine is skipped */
{
base = lyapzero + vec_dim*i;
y_temp[base]
= -2*y[0]*y[ base ]
+ C1*y[ base + 1];
y_temp[base+1]
= y[ base ];
}
for (i = 0; i < lyapzero + vec_dim*num_lyap; i++)
y[i] = y_temp[i] ;
}
ADDING A DIFFERENTIAL EQUATION
The program knows that a differential equation is going to be
used when it finds the differential equation step size has been set,
that is the variable "step",e.g.: step = 1./100.;
For differential equations, the program must be told in the
initialization file what the name of the vector field is, DELorenz in
this case. Thus the function pointer DE is set in the Lorenz
initialization with the line:
DE = DELorenz;
The pointer DE is called by rungekutta(), or any other solver used.
The Lorenz vector field is defined as follows; Y[] denotes the
position at which the vector field will be evaluated and k[] is the
vector field. Notice that there are 4 coordinates, one of which is time.
Just as with y[], k[] must have coordinates for phase space, time, and
all lyapunov vectors.
DELorenz(k,Y)
double Y[],k[];
{
int i,base;
double x1,y1,z1;
double J00,J01,J10,J11,J12,J20,J21,J22;
x1 = Y[0];/* this and the following 2 lines are for speeding
computation ever so slightly and for increasing
readability */
y1 = Y[1];
z1 = Y[2];
k[0] = -sigma*(x1 -y1);
k[1] = rho*x1 -y1 -x1*z1;
k[2] = x1*y1-beta*z1;
k[3] = 1; /* time */
if (num_lyap > 0) /* the following are the partial derivatives
of the vector field */
{
J00 = -sigma;
J01 = sigma;
J10 = rho-z1;
J11 = -1.;
J12 = -x1;
J20 = y1;
J21 = x1;
J22 = -beta;
}
for (i = 0; i < num_lyap; i++)
{
base = lyapzero + vec_dim*i;
k[base + 0] =
J00*Y[base] + J01*Y[base+1];
k[base + 1] =
J10*Y[base] + J11*Y[base+1] + J12*Y[base+2];
k[base + 2] =
J20*Y[base] + J21*Y[base+1] + J22*Y[base+2];
}
}
initLorenz()
{
int DELorenz();/* this is not needed if DELorenz() is defined
in the same file as this routine and is earlier in the
file than this routine. */
vec_dim = 3; /*the dimension of the
Lyapunov vectors = phase space dim*/
num_lyap = 0; /*the number of Lyapunov numbers
to be computed <= vec_dim */
lyapzero = 4;/* y[lyapzero] is the zeroth coord of the
zeroth lyapunov vector */
X_upper = 25.0; /* x scale */
X_lower = -25.;
Y_upper = 60.; /* y scale */
Y_lower = 0.;
Y_coord = 2;
sigma = 10;
rho = 28.0;
beta = 8./3.;
step = 1./100.;
DE = DELorenz; /* DE is a pointer to a function */
}